pacman::p_load(tidyverse,data.table,broom,openxlsx,latex2exp,ggpubr)
rm(list=ls())
`%!in%` = Negate(`%in%`)
################################################################################

###Spatial units
df_sa_units = setDT(read.xlsx('aux_params/parameters.xlsx', sheet='spatial_units'))

###Emissions intensities
df = setDT(read.xlsx('aux_params/parameters.xlsx',  sheet = "emissions_nluc"))
setnames(df, old = c('beef_nluc_per_t', 'crops_nluc_per_t'),  new = c('beef_E_t','crops_E_t'))
df = merge(df[, list(spatial_id_name,beef_E_t,crops_E_t)], df_sa_units[,list(spatial_id_name, region, frontier)], by='spatial_id_name', sort=F)
df_nluc = df
df = setDT(read.xlsx('aux_params/parameters.xlsx',  sheet = "emissions_luc"))
setnames(df, old = c('tCO2_per_ha_stock'),  new = c('E_ha'))
df = merge(df[, list(spatial_id_name,E_ha)], df_sa_units[,list(spatial_id_name, region, frontier)], by='spatial_id_name', sort=F)
df_luc = df
df=merge(df_nluc,df_luc, by=c('spatial_id_name','region','frontier'), sort=F)
df_emissions=df[,list(spatial_id_name, beef_E_t, crops_E_t, E_ha, frontier)]

for(counterfactual in c('dt','tariff') ){
    
    ##### Changes
  
    #LF
    df = read.csv(file = paste0('results_ic/lf/Q_ij.csv'), header=FALSE)
    setnames(df, old=c('V1','V2','V3','V4'), new=c('Domestic','EU','Asia','ROW') )
    df_lf_Q = df
    #CF
    df = read.csv(file = paste0('results_ic/',counterfactual,'/Q_ij.csv'), header=FALSE)
    setnames(df, old=c('V1','V2','V3','V4'), new=c('Domestic','EU','Asia','ROW') )
    df_cf_Q = df
    #CF-LF
    dQ = df_cf_Q-df_lf_Q
    dQ_b = dQ[1:(nrow(dQ)/2),]
    dQ_c = dQ[(1+(nrow(dQ)/2)):nrow(dQ),]
    
    #LF
    df = setDT(read.csv(file = paste0('results_ic/lf/L.csv'), header=FALSE))[,list(V1,V2)]
    setnames(df, old=c('V1','V2'), new=c('beef','crops') )
    df$agriculture = df$beef + df$crops
    df_lf_L = df
    #CF
    df = setDT(read.csv(file = paste0('results_ic/',counterfactual,'/L.csv'), header=FALSE))[,list(V1,V2)]
    setnames(df, old=c('V1','V2'), new=c('beef','crops') )
    df$agriculture = df$beef + df$crops
    df_cf_L = df
    #CF-LF
    dL = setDT(df_cf_L-df_lf_L)
    
    ##### Realized yields
    
    dfb = df_lf_Q[1:(nrow(df_lf_Q)/2),]
    A_b = (dfb$Domestic+dfb$EU+dfb$Asia+dfb$ROW)/df_lf_L$beef
    dfc = df_lf_Q[(1+(nrow(df_lf_Q)/2)):nrow(df_lf_Q),]
    A_c = (dfc$Domestic+dfc$EU+dfc$Asia+dfc$ROW)/df_lf_L$crops

    ##### Emissions
    
    #NLUC
    dE_b_nluc = dQ_b*df_emissions$beef_E_t
    dE_c_nluc = dQ_c*df_emissions$crops_E_t
    
    #LUC
    df = setDT(dQ_b/A_b) + setDT(dQ_c/A_c)
    df = df[, c('TOTAL'):=list(Domestic+EU+Asia+ROW)][, c('share_Domestic','share_EU','share_Asia','share_ROW'):=list(Domestic/TOTAL, EU/TOTAL, Asia/TOTAL, ROW/TOTAL)]
    df = df[, list(share_Domestic,share_EU,share_Asia,share_ROW)]
    df = df*dL$agriculture
    df = df[, c('TOTAL'):=list(share_Domestic+share_EU+share_Asia+share_ROW)]
    df=df*df_emissions$E_ha
    setnames(df, old=c('share_Domestic','share_EU','share_Asia','share_ROW'), new=c('Domestic','EU','Asia','ROW'))
    dE_luc=df[, list(Domestic,EU,Asia,ROW)]
    
    #Total
    df = setDF(dE_b_nluc) + setDF(dE_c_nluc) + setDF(dE_luc)
    df = setDT(df)[, lapply(.SD, sum, na.rm=TRUE)]
    df$nonEU = df$Domestic + df$Asia + df$ROW
    df$TOTAL = df$Domestic + df$EU + df$Asia + df$ROW
    df = setDT(gather(df, destination, delta, Domestic:TOTAL, factor_key=TRUE))
    
    if(counterfactual=='tariff'){plot_title='A. Incomplete regulation: unilateral EU tax'}
    if(counterfactual=='dt'){plot_title='B. Complete regulation: domestic downstream tax'}
    
    figure = ggplot(df[destination %!in% c('nonEU')], aes(x=destination, y=delta/1000000, fill = ifelse(destination == "TOTAL", "Highlighted", "Normal"))) +
      geom_bar(stat='identity', position="dodge", color='black', linewidth=0.1, alpha=1, width=0.5) +
      labs(x=TeX('destination market'), y=TeX('$\\Delta$ Mt $CO_2$e'), title=plot_title)+
      geom_hline(yintercept=0, color='black', linewidth=0.25)+
      theme_minimal()+
      theme(aspect.ratio = 0.6,
          plot.title = element_text(size=12, vjust=3),
          axis.title.y = element_text(size=9, hjust=0.9),
          axis.title.x = element_text(size=9, margin = margin(t = 10)),
          axis.text.x = element_text(size=9),
          axis.text.y = element_text(size=7),
          panel.grid.minor = element_blank(),
          legend.title = element_blank(),
          legend.text = element_text(size=6),
          legend.key.height=unit(0.3,"cm"),
          legend.position = 'none')+
    scale_fill_manual(values = c("Highlighted" = "gray25", "Normal" = "lightgray"))

    if(counterfactual=='tariff'){figure_tariff=figure}
    if(counterfactual=='dt'){figure_dt=figure}
    
}

figure=ggarrange(figure_tariff,figure_dt,nrow=1,ncol=2)
ggsave(paste0('../../output/figures/appendix/figure7_leakage.pdf'), width = 9, height = 3.5)
ggsave(paste0('../../output/figures/appendix/figure7_leakage.eps'), width = 9, height = 3.5)

